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We present the complete set of symmetric and antisymmetric {edge and corner) surface modes in 
finite one- and two-dimensional arrays of waveguides. We provide classification of the modes based 
on the anti-continuum limit, study their stability and bifurcations, and discuss relation between sur- 
face and bulk modes. We put forward existence of surface breathers, which represent two-frequency 
modes localized about the array edges. 

PACS numbers: 42.65.Tg, 42.65.Sf, 42.65.Wi 



I. INTRODUCTION 



II. SURFACE AND BULK MODES 



A. The model and terminology 
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Waves at surfaces and interfaces are known to exhibit 
peculiar properties. Localized electronic states at a crys- 
tal edge, discovered by Tamm were the first example 
of such phenomena. Later on it was found, that surfaces 
and interfaces are able to sustain localized waves which 
attracted a great deal of attention in different areas of 
physics, particularly due to variety of practical appli- 
cations, like plasmonic waveguides 0, sensors Q, etc. 
Recently it was predicted theoretically 'i'l and observed 
experimentally ,5] , that at the edge of a semi-infinite one- 
dimensional (ID) array of nonlinear waveguides there can 
exist discrete surface solitons. Modes localized at finite 
distances from the edge were considered in ^ and surface 
gap solitons between uniform media and periodic lattice 
were reported in Q • Different type of surface modes in 
2D arrays were studied in [1, Q . 



Structures with two surfaces give rise to novel proper- 
ties of surface modes. For example, interaction of the two 
surface polaritons, supported by each surface of a metal- 
lic film, results in creation of symmetric and antisymmet- 
ric modes [lo| . which in their turn originate polariton- 
assisted extraordinary transmittancy of the film [ llj . 



In this paper we describe discrete (edge and corner) 
surface modes in finite ID and 2D arrays of nonlinear 
waveguides. We show that they can be classified on 
the basis of the anti-continuum (AC) limit, similarly to 
the classification of intrinsic localized modes introduced 
in p^ . and in this way the complete families of modes 
can be identified. We show that surface modes can bi- 
furcate either with other surface or with bulk modes and 
study the mode stability. We also report a new type of 
the surface excitations - surface breathers - which repre- 
sent two-frequency excitations localized in the vicinity of 
the array edges. 



We start with a finite array of M nonlinear waveguides 
described by the discrete nonhnear Schrodinger (DNLS) 
equation [Tj] 



M 



iqn + y ^ ((^n'.n+l + Sn' ,n-l)qn' + O" Qn = 0. (1) 



ri' = l 



Here g„ = dqn/dC, ^ is the propagation coordinate, q„ 
is the dimensionless field amplitude in the n-th waveg- 
uide {n = 1,...,M), a — 1 and a = —1 stand for fo- 
cusing and defocusing nonlinearities. We concentrate on 
modes having definite parity imposing (/„ — QM+i-n for 
symmetric and Qn = —QM+i^n for antisymmetric modes. 
Eq. II]) possess two integrals of motion: the Hamiltonian 

CT V— ^Af 



and the 



total power P = Yln=i \lnf- 

It worth to emphasize that the DNLS equation ([T]) is a 
widely used model in the condensed matter physics [3| 
and in the theory of Bose-Einstein condensates loaded 
in optical lattices [l^, what makes the results reported 
below to be relevant for a rather wide class of the phe- 
nomena of the nonlinear physics of periodic and discrete 
structures. 

Like in the well studied infinite case [IB], station- 
ary modes of Eq. ([1]) are searched in the form (?„(C) = 
Qn exp(— jAC), where A is the propagation constant, and 
the resulting equations are 



M 

E 

ji'=i 



)Qn' + (tQI - 0. (2) 



The solutions for a — ±1 are connected by the following 
symmetry reduction 12]: if the Qn is a solution of ([!]) 
for a definite A and a = +1, then (— 1)"(5„ is a solution 
for —A and <t — —1. 

In order to describe the whole diversity of solutions 
and to classify them [13], one has to consider the AC 
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limit [T3| ■ To this end we rewrite Eq. ([T]) in terms of the 
rescaled stationary amplitudes Vn = Qn/ ^ 

M 

<^\M'>Jn{vl + CTS) + X! ("^"'.n+l + Sn',n-l)Vn' = 0, (3) 
n' = l 

where s = sign(A), and consider the limit |A| — > oo. 
In this limit t;„ become independent and for the case 
as = —1 acquire one of the three values: u„ = — 1, 
Vn — 0, and w„ = +1 (1 < n < AI). Thus, in the 
AC hmit there exists N, = (3[(m+i)/2] _ symmet- 
ric and Na = (s'*^/^' — l) /2 antisymmetric modes (here 
square brackets signify the integer part), and each mode 
can be coded [l^ by a sequence of M symbols — , 0, and 
+. The coding corresponds to the limit of the infinite 
power: i.e., in particular, "0" does not refer to the zero 
intensity of a waveguide at a finite input intensity (see 
Fig. [21 and discussion below). As an example, an array of 
three waveguides has four symmetric: {0 + 0}, {+ + +}, 

{H h}, {+0+} modes, and one antisymmetric mode 

{+0— }. A sequence, consisting of symbols {+,0,—}, is 
termed a " word' , a word having only zeros is referred 
to as "empty", and a word having no (for example 

{H h}) is called a "simple" word. Number of symbols 

in a word is called the length of the word. 

The first important property of the coding steams 
from the analytical continuation of the AC limit to 
|A| > Aac, where Aac is a constant (in the case of the 
infinite array Aac ~ 5.4533 [l^). This means that all the 
words exhaust all possible modes of a finite array exist- 
ing for A > A* . Second, AC limit allows one to introduce 
a definition for a surface mode as a code consisting of 
two simple words separated by an empty word, the lat- 
ter having the length not less than the lengths of each of 

the simple words. For example {H hOOOO -\ h} is 

a surface mode of an array of 10 waveguides. All other 
modes will be referred to as bulk modes. The introduced 
terminology, being mathematically well defined, has rel- 
ative physical meaning for finite A: surface modes can 
bifurcate with bulk modes, the both acquiring identical 
shapes in the bifurcation point. 

B. Bulk modes 

In the linear limit P ^ (or formally a 

0) the Eq.® possesses M eigenvalues Aq™^ = 
—2 cos [ttto/ {M + 1)] corresponding to eigenmodes 

«=sin(^), .n.l,...,M. (4) 

Thus one can expect that M bulk modes have linear 
limit and thus do not possess an intensity threshold of 
excitation: for these modes, when A — > Aq™'*, the power 
p(™) of m-th mode (m = 1, . . . , M) approaches zero (see 
FigH])- In order to determine the dependence p('")(A) 




FIG. 1: Dependence of power P vs the propagation constant 
A for the bulk symmetric and antisymmetric modes in the case 
(7 = 1, M = 6, which do not have the excitation threshold: 1 - 

mode {00-h-hOO}, 2 - mode {0-hOO-O}, 3 - mode {+0 0-f }, 

4 - mode {+0 - -1-0-}, 5 - mode {-I h H h}, 6 - mode 

{H 1 1 — }. The bifurcation points with asymmetric modes 

are depicted by filled circles. 

near the linear limit we follow the standard perturbation 
technique (see, e.g., [l^l), and look for a solution of ^ 
in a form of the series 

Qn^eQ^^!+e'Q^^J+o{e% (5) 

A = A{,")+e2A(")+o(e2), (6) 

where we have introduced the small parameter e = 
^/2P/{M + 1) < 1. Substituting the above expansions 
into Eq. ^ and gathering the terms of the same order in 
e, we rewrite Eq.(l2]) in the form of a set of equations: 

M 

A^^QS + E + ^n,n-^)Q'^ = F^. (7) 
n' = l 

Here F^^) = 0, = -A^^Q^™) - a (q("J)'. As 

it is clear Eq.Q for j = describes a linear eigenmode 
and therefore is automatically satisfied, while consider- 
ing solvability conditions for j = 2 (it is equivalent to 
orthogonality ^^2'^'' ^^d Qg™"*, we obtain the corrections 
to the eigenvalues written in the form 

It follows from Eq.® that each of M linear modes 
possesses its unique small-amplitude nonlinear analogue 
(from all diversity of nonlinear modes only given M 
bulk modes have small-amplitude solution). This also 
proves that no linear surface mode exists (what corrob- 
orates with the earlier findings for a semi-infinite ar- 
ray Moreover, from Eq.® it follows, that in the 
small-amplitude limit these modes are characterized by 
the linear dependence of the mode total power upon the 
propagation constant: 

=a{2M + 2) -2 . (9) 

3 + 0,n,{M+l)/2 

From Eq.® it is clear, that powers of these modes are 
decaying for cr = 1 and increasing for a = —I functions 
of A. To single out the obtained modes in what follows 
they are referred to as quasi-linear modes. 
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C. Bifurcations of quasi-linear modes 

When the power increases one of two scenario of mode 
transformations is possible: either the branch of the 
quasi-hnear mode smoothly tends to a uniquely defined 
AC limit or at some A* (alternatively P*) it bifurcates 
giving origin to some new solutions. As it is clear, each 
mode can bifurcate either with a mode of the same sym- 
metry or with an asymmetric mode (which in principle 
does not possess any symmetry, but in the bifurcation 
point acquires the given symmetry). If the former even 
takes place then the quasi-linear mode is naturally classi- 
fied by its AC limit. While analytical description of other 
cases we leave for further studies, we mention that all nu- 
merical simulations we performed with a finite number of 
waveguides have shown that no bifurcations of the quasi- 
linear modes with modes of the same symmetry occurs: 
only asymmetric modes bifurcate from the quasi-linear 
ones. This allows us to use for the latter modes the classi- 
fication determined by their symmetric AC limit (i.e. by 
a word describing symmetric ramification of the mode). 
Namely this notations are used in the figure captions and 
in the text whenever we speak about quasi-linear modes. 

To determine numerically the bifurcation points of the 
quasi-linear modes we consider the continuation of the 
mode by the parameter A (notice that contrary the stan- 
dard approach [l^, now we " move" along the branch 
outwards the linear Hmit) 



dA 

where Q — col {Q,. 
F(Q,A) - 



[DqF] 



9F(Q,A) 
dX 



(10) 



/ M \ 

col XQn + ^ ((5„',„+l -I- Sn',n-l)Qn' + CtQI , 

and the entries of the three-diagonal matrix DqF are 



(DqF) 



(A -f- 3aQl) 



-'n,n' + l ■ 



(11) 



This continuation of the quasi-linear modes is possible, 
while the matrix -DqF is invertible, i.e. its determinant 
is not equal to zero. Thus the bifurcation point is deter- 
mined by the equation D = D{X) = DetlDgFj = 0. 

As we already mentioned M quasi-linear modes can 
bifurcate with asymmetric modes what is illustrated in 
FigH] obtained for M = 6. There all symmetric modes 
possess additional bifurcation points, denoted by filled 
circles. Although the consideration of modes without 
symmetry is beyond the scope of this paper, we now an- 
alyze analytically the simplest case of M = 2 and a = 1 
allowing the trivial solution. In that case there exists one 
antisymmetric mode Qi 2 = 1 ^ A (with code {-I — }) and 
one symmetric mode Ql 2 = — 1 — A (with code {++}). 
Substituting the expressions for the field distribution of 
the antisymmetric mode {H — } into (jlip . we obtain, that 



D{X) = only when A = Aq = 1 (where the mode is 
born). For the symmetric mode {H — h} one verifies that 
similar to the previous case D{X) = at the point of lin- 
ear limit A = Aq^^ = —1, but also at A = A, = —2, where 
the mode {++} bifurcates with mode {+0}. The mode 
{+0} (for which 2 = -A/2 ± x/X^/i- 1) exists only 
when A < A*, and at the point A* its has the same field 
distribution as mode {++} (so, we have a pitchfork-type 
bifurcation) . 




FIG. 2: P vs X for fundamental (a), (d), and twisted (g) 
modes. In panels b,c,e,f examples of the symmetric (A,B,E,F) 
and antisymmetric (C,D,G,H) fundamental surface modes, as 
well as their classifications, are shown for M = 7 (panels b,c) 
and A/ = 6 (panels e,f). In panels h,i symmetric (A,B) and 
antisymmetric (C,D) twisted surface modes for M = 7 are 
shown. In all examples (7 = 1. For comparison, in panels a,d, 
and g the dashed lines represent doubled powers of the same 
modes in a semi-infinite waveguide array. 



D. Surface modes in one-dimensional finite lattice 

Turning now to the analysis of surface modes, by anal- 
ogy with an infinite array [Toj . one can distinguish fun- 
damental modes having in-phase distribution of the field 
near the edges and twisted surface modes, having out- 
of-phase fields in the two waveguides bordering an edge. 
As an examples, in Fig. [2] we show the mode patterns 
for arrays of M = 6 and M = 7 waveguides. All the 
modes shown require a threshold power to be excited. 
Two different symmetric and antisymmetric fundamen- 
tal surface modes bifurcate with each other at aI^''*'' 
and aI^^"-* (the codes of such modes are indicated in 
the respective panels of Fig. [5]). The pairs of symmet- 



ric modes, (A,B) and (E,F), bifurcate in A: 



(7s) 



-3.005 

and aI^^'' ~ —3.115, while the bifurcation points of 
antisymmetric modes (C,D) and (G,H), are given by 

2.67. We observe that 



a1'"' « -2.83 and aI'^^ 



aI*^'*'' < aI*^'^'' and thus the antisymmetric modes are 
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excited at lower field intensities. Comparing the bifur- 
cation points of the modes with different M we also ob- 
serve that for large enough arrays, i.e. at Af ^ oo, the 
modes are transformed in the conventional surface modes 
of a semi-infinite array. In this limit distinction between 
symmetric and antisymmetric modes disappears. Thus, 
presence of two boundaries of a finite array essentially 
modifies a surface mode. The physical reason for this 
that in a vicinity of the bifurcation point the modes are 
weakly localized near the array edges, and the field in 
waveguides near the center of the array is non negligible 
(see Fig 12). This leads to interaction between the modes, 
supported by the two edges, what in its turn modifies the 
patterns. At large values of |A| the modes are strongly 
localized near the edges and interaction between them is 
weak, what results in the identical asymptotic behavior 
of symmetric and antisymmetric modes in the AC limit 
clearly seen in Fig. [5] a, d. 



The code of the twisted mode includes more nonzero 
symbols in comparison with the code of the fundamen- 
tal mode. So, since symbol "0" in the code of certain 
mode signifies zero field amplitude in the correspondent 
waveguide in the AC limit, the power of twisted mode in 
AC limit should be higher than the power of fundamen- 
tal mode. Nevertheless as shown in Fig.[^g, even in the 
vicinity of the bifurcation point twisted modes are excited 
at higher intensities than the fundamental modes. Also 
due to higher field in the center of the waveguide array (in 
comparison with the fundamental mode case) the inter- 
action between two edges of the array is stronger, what is 



(7a)_^(7s) 



> A 



(7a) 
[{7a) 



-A1'^) (the 



-3.665 and A!,;"^ « -2.83). 



expressed by the relation A 
bifurcation points are A^^''"' 
The main peculiarity of the diagram of the twisted modes 
is that surface modes A and C bifurcate with the bulk 
modes B and D. 



To study linear stability of the modes, we follow the 
standard steps analyzing the eigenvalue problem lin- 
earized about the mode, where /3 is the spectral parame- 
ter such that Im(/3) < corresponds to a linearly unsta- 
ble mode. The imaginary parts of P for the fundamental 
modes B, D, F, and H decrease with A —oo, i.e. the 
modes are unconditionally unstable (similar behavior was 
reported in [1^). Symmetric fundamental modes A and 
E are unstable, but Im(/3) increase with the decreasing of 
A. The antisymmetric fundamental modes C and G are 
unstable in the vicinity of the bifurcation point, but are 
stabilized at A < -3.13 (for M = 7) and A < -3.22 (for 
M — 6). Stability analysis of the twisted modes shows, 
that symmetric modes A, B, and the antisymmetric mode 
D are unconditionally linearly unstable, while the anti- 
symmetric mode C is stable at A < —3.87 (these results 
corroborate with the results of [2l[ for the discrete modes 
of infinite array). 



III. SURFACE BREATHERS 

Now we consider a novel type of localized modes - 
surface breathers. Antisymmetric surface breathers can 
be constructed analytically for an array of M = 5 sites. 
In this case, Eq. ([T]) possesses a solution of the following 
type (found using the known dynamics of a dimer (2^ ) 



cn(C(C-Co)/^,fc), < A: < 1 
dn(C(C-Co),l/fc), k>l 



(12) 



where cn and dn are Jacobi elliptic functions, z = 
2(|9iP — \q2\'^)/P is the intensity contrast, k ~ 
C/{V^g{P)) is the elliptic modulus, 

Co = fcF(arccos(z(0)F/(4C)),fc)/C 

F{(j), k) is the incomplete elliptic integral of the first kind, 

= - H/2- P^/16-2, 

and g = + 2H + 4)1/-^. When fc > 1, the function 

z(C) « (4C/P) [1 - sin2(C(C - Co))/2fc2] ^ 

describes a mode concentrated near the edges n = 1, 5 
(see Fig. [3]) and oscillating with the period 2K{l/k)/C 
{K{k) is the complete elliptic integral of the first kind). 
In the limit fc — > cxd (taken at a constant power P) the 
Hamiltonian achieves it minimal value H — — P^/4 — 2 
and the intensity contrast becomes a constant. In that 
case the surface breather transforms into a fundamen- 
tal surface mode (e.g. the breathers in Fig. [3] a and b 
transform into the surface modes with P ~ 7, X — —3.5, 
H = -14.25 and P = 12, A = -6, iJ = -38, respec- 
tively). Thus, a surface breather can be excited from 
the fundamental surface mode by small detuning of its 
Hamiltonian from the minimum. We used this idea to 
construct the antisymmetric surface breathers for arrays 
of M = 7 and M = 8 waveguides (Fig.[3]c and d), which 
can not be constructed analytically. Notice that when 
the Hamiltonian possesses its minimal value, the surface 
breathers of Fig. [3] c,d transform into stationary modes 
with M = 7, X = -6,H = -38.185 and M = 8, A = -6, 
H = —38.192, respectively. 

We tested the stability of breathers by direct numeri- 
cal solution of differential Eq. ([T|), perturbing the initial 
profiles by noise with an amplitude of order of 10% of 
Qn in each waveguide. The breather, depicted in Fig. [3] 
a has shown unstable behaviour, while breathers of Fig. [3] 
b-d demonstrated a stable one. Like in the case of the 
surface modes increasing of |A| results in stabilization of 
a surface mode. 

The case M — A also allows for analytical construc- 
tion of symmetric (b ~ 1) and antisymmetric (b = 
-1) breathers. Now z = zi + 6/'(zi)/[24p(P(C - 
Co)/4;52,53) - f"{zi)], where zi is a root of the poly- 
nomial f{z) — a4 + a^z + a2Z^ + aiz^ — z^, 04 — 
16{4:~4:H^ /P^-H-P^/ 16) /P'^, 03 = 6ib{H+P'^/8)/P^, 
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FIG. 3: (Color online) Surface breather for (a) M = 5, P = 7, 
H = -14.175; (b) M = 5, P = 12, H = -37.602; (c) M = 7, 
P = 12.021, H = -37.463 and (d) M = 8, P = 12.022, 
H — —37.471, excited by noise. 



02 = -8(10+2i/+P74)/P^ ai = 8b/ P, and p(a;; 52, 53) 
is the Weierstrass elliptic function with 52 = ^l/l^ ~ 
0103/4 — 04 and 53 = 010203/48 + 0!- O2/2I6— (01/I6 + 
02/6)04. The character of oscillations of the field along 
the waveguides is determined by A = (/f- 27(7|: for A ^ 
the solutions are oscillatory about a nonzero average, but 
for A = the solutions are aperiodic. 



IV. CORNER AND EDGE MODES 

Now we discuss surface modes of a 2D finite M x M 
array, where each waveguide is coupled with the nearest 
neighbors. The system is described by coupled 2D DNLS 
equation 



M 



m' — l 

M 

{^n',n+l + Sn',n-l) ^ri',™ + CT \qn,m\ qn,m — 0.(13) 



Now the total power is given by P = Y,n,m=i 19", ml 
and the symmetry reductions of the stationary modes 
9n,m(C) = Qn,mexp(-'iAC) are as follows: if the Qn,m 
is a solution of p3)) for a definite A and ct = +1, then 



(— l)"+'"Q„^m is a solution for —A and a = —1. 

Similar to the ID case, each mode can be coded on a 
2D map by an array of M x M symbols — , 0, and + (see 
examples in Fig.|4|), corresponding to its AC limit {P — > 
oo). Splitting the array into blocks, we call "empty" 
a block consisting of all zeros and " simple" a simply- 
connected block having no zeros. Now we define a corner 
mode as a mode on a square array consisting of only sim- 
ple blocks at the corners, separated by empty blocks of 
higher dimensions (see an example in FigH^). Being in- 
terested only in modes of a definite symmetry we impose 
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FIG. 4: Examples of the fully symmetric corner (a) and edge 
(b) modes. The dashed lines outline empty and simple blocks. 



additional constrains: qn,m = QM+i-n.m = qnM+i-m 
(for fully symmetric modes), g„,™ = -qM+i-n,r,i = 
qn.M+i-m (for symmetric-antisymmetric modes) and 
<?ri,m = -(jM+i-n,™ = -Q„,M-i-i-m (for fuUy antisym- 
metric modes). Similarly we identify an edge mode, 
whose code consists of simple blocks bordering edges of 
the array, but separated from each other and from the 
corners by empty blocks (see the example in Fig. HJd). 
A symmetry constrains for the edge mode are as fol- 
lows: g„,,„ = g„i,„ = qM+i-m,M+i-n (for Mly sym- 
metric modes), qn,yn = q-m.n = -qM+i-m,M+i-n (for 
symmetric-antisymmetric modes) and qn,m = —qm,n = 
~qM+i-m.M+i-n (for fuUy antisymmetric modes). 

We restrict the consideration to the lowest-power fun- 
damental modes, i.e. modes coded by one-site simple 
blocks. We found that similar to the ID case both cor- 
ner and edge fully antisymmetric modes (C and F in 
Fig. [5l correspondingly) require lower threshold power 
of excitations that other types of the modes, while the 
fully symmetric modes (A and D in Fig. [51 correspond- 
ingly) are excited at higher powers (see upper panel in 
Figin]). At the same time the distinction between the 
properties of the fully symmetric, semi-symmetric and 
fully antisymmetric edge modes is stronger than the dis- 
tinction between the properties of similar types of corner 
modes. This occurs due to the smaller distance (and, 
hence, stronger interaction) between the excitation cen- 
tered at four edges in comparison with excitations cen- 
tered at four corners (compare, e.g. field patterns for 
modes A and D, B and E, C and F in FigO. The stabil- 
ity analysis shows that fully symmetric and symmetric- 
antisymmetric modes A,B,D,E are unstable, but their 
instability increments Im(/3) < increase as A decreases. 
The fully antisymmetric modes C and F are stabilized at 
A < —4.74 and A < —5.4, respectively. 



V. CONCLUSION 

We have reported the complete families of (edge and 
corner) surface modes in arrays of ID and 2D waveguides, 
presenting the classification exhausting all possible sta- 
tionary excitations. It has been shown that the surface 
modes belong to one-parametric branches of solutions, 
which bifurcate either with other surface or with bulk 
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FIG. 5: (Color online) The total intensity P vs the propa- 
gation constant A (upper panel) for corner (solid lines) and 
edge (dashed lines) modes. Examples of the patterns of the 
corner (panels A-C) and edge (panels D-F) modes are shown 
for M — 7 and cr = 1. A and D are fully symmetric modes, B 
and E are symmetric-antisymmetric modes, and C and F are 
fully antisymmetric modes. 



modes, when the total power is properly changed. We 
have also found two-periodic modes, surface breathers, 
whose intensity periodically oscillates staying localized 
about the edge of the array, and described a way of ex- 
citations of breathers starting with the respective sur- 
face modes. The reported solutions, being well local- 
ized at the two (or four in the 2D case) surfaces and 
at the same time showing nonzero intensity in the bulk 
of the array, can be of significant practical importance, 
by analogy with the polaritons assisting extraordinary 
transmittancy . Meantime we emphasize a number of 
open questions which were left unanswered by the present 
research. Among those we mention thorough study of 
the linear stability of the surface breathers, mathemati- 
cal justification of the complete classification of the two- 
dimensional modes starting with the AC limit. These 
issue as well as specific practical outputs will be address 
elsewhere. 



Acknowledgments 

YVB. was supported by the FCT Grant No. 
SFRH/PD/20292/2004. VVK. acknowledges support of 
the Sccrctaria de Stado de Universidades e Investigacion 
(Spain) under Grant No. SAB2005-0195. The work was 
supported by the FCT and European program FEDER 
(Grant No. POCI/FIS/ 56237/2004). 



[1] 
[2] 

[3] 

[4] 
[5] 
[6] 

[7] 

[8] 



[9] 

[10] 

[11] 
[12] 



I. E. Tamm, Z. Phys. 76, 849 (1932). 

P. Berini, R. Charbonneau, N. Lahoud, and G. Mattiussi, 

J. Appl. Phys. 98, 043109 (2005). 

see e.g. A. G. Brolo, R. Gordon, B. Leathem, and K. L. 

Kavanagh, Langmuir 20, 4813 (2004). 

K. G. Makris, et al Opt. Lett. 30, 2466 (2005). 

S. Suntsov, et al Phys. Rev. Lett. 96, 063901 (2006). 

M. L Molina, R. A. Vicencio, and Y. S. Kivshar, Opt. 

Lett. 31, 1693 (2006). 

Y. V. Kartashov, V. A. Vysloukh, and L. Torner, Phys. 
Rev. lett. 96, 073901 (2006). 

K.G. Makris, J. Hudock, D.N. Christodoulides, G.L 
Stegeman, O. Manela, and M. Segev, Opt. Lett. 31, 
2774 (2006); H. Susanto, P.G. Kevrekidis, B.A. Mal- 
omed, R. Carretero-Gonzalez, and D.J. Franzeskakis, e- 



print nlin.PS/0607063 



R. A. Vicencio, S. Flach, M. I. Molina, and Y. S. Kivshar, 
e-print cond-mat/0610049 

R. E. Camley and D. L. Mills, Phys. Rev. B 29, 1695 
(1984). 

T. W. Ebbesen, et Nature 391, 667 (1998). 

G.L. Alfimov, V.A. Brazhnyi, and V.V. Konotop, Phys- 

ica D 194, 127 (2004). 



[13] D. N. Christodoulides and R. I. Joseph, Opt. Lett. 13, 
794 (1988). 

[14] A. Scott, "Nonlinear science: emergence of coherent 
structures." (Oxford University Press, 1999) 

[15] V. A. Brazhnyi and V. V. Konotop, Mod. Phys. Lett. B 
18, 627 (2004). 

[16] see e.g. D. Hennig, and G. Tsironis, Phys. Rep. 307, 333 

(1999); P. G. Kevrekidis, K. 0. Rasmussen, and A. R. 

Bishop, Int. J. Mod. Phys. B15, 2833 (2001). 
[17] R. S. MacKay and S. Aubry, Nonlinearity 7, 1623 (1994). 
[18] R. Bellman, "Introduction to Matrix Analysis." 

(McGraw-Hill Education, 1970) 
[19] S. Darmanyan, A. Kobyakov, and F. Lederer, JETP 86, 

682 (1998) 

[20] Y. V. Kartashov, V. A. Vysloukh, D. Mihalache, and 

L. Torner, Opt. Lett. 31, 2329 (2006). 
[21] D.E. Pelinovsky, P.G. Kevrekidis, D.J. Frantzeskakis, 

Physica D 212, 1 (2005) 
[22] V. M. Kenkre and D. K. Campbell, Phys. Rev. B 34, 

4959 (1986); S. Raghavan, A. Smerzi, S. Fantoni, and 

S. R. Shenoy, Phys. Rev. A 59, 620 (1999). 



